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Abstract 


In this article, we constructed a numerical scheme for singularly perturbed 
time-delay reaction-diffusion problems. For the discretization of the time 
derivative, we used the Crank-Nicolson method and a hybrid scheme, which 
is a combination of the fourth-order compact difference scheme and the cen- 
tral difference scheme on a special type of Shishkin mesh in the spatial di- 
rection. The proposed scheme is shown to be second-order accurate in time 
and fourth-order accurate with a logarithmic factor in space. The uniform 
convergence of the proposed scheme is discussed. Numerical investigations 
are carried out to demonstrate the efficacy and uniform convergence of the 
proposed scheme, and the obtained numerical results reveal the better per- 
formance of the present scheme, as compared with some existing methods 
in the literature. 
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1 Introduction 


Delay differential equations play a crucial role in the mathematical modeling 
of various practical phenomena and are widely applicable in fields such as 
biosciences, control theory, economics, material science, medicine, robotics, 
etc. [17, 16, 15]. If we restrict the delay differential equations to a class 
in which the highest order term is multiplied by a small parameter ¢, then 
we call it singularly perturbed delay differential equations of the retarded 
type. Nowadays, there has been a growing interest in the study of singularly 
perturbed delay differential equations due to their occurrence in many prac- 
tical models, for instance, Hutchinson’s equation which is a model for the 
evolution of the population in mathematical ecology|1] and, the Wazewska- 
Czyzeska and Lasota equation that describes the survival of red blood cells 
in animals [27]. 

The solutions of singularly perturbed time-delayed problems are funda- 
mentally different from those problems without time delay because small 
lags can have large effects. The solution of delay differential equations re- 
quires knowledge of both the current state and the state at a certain time 
in the past. Finding the solution to singularly perturbed delay differential 
equations using classical numerical methods fails to give stable and accurate 
results because of the presence of the perturbation parameter e. 

Many researchers have treated time-dependent singularly perturbed parabolic 
partial differential equations,with or without time-delay: to mention some 
[3, 4, 14, 12, 2, 10, 20, 21, 22, 23, 28]. However, time-delayed reaction diffu- 
sion problems have not been extensively investigated. Authors in [13] solved 
singularly perturbed time-delay parabolic partial differential equations of the 
reaction-diffusion type. The problem is discretized by a hybrid scheme on a 
generalized Shishkin mesh in the spatial direction and the implicit Euler 
scheme on a uniform mesh in the time direction, and to increase the order 
of convergence in the time direction, Richardson extrapolation is applied. 
The authors in [12] studied a similar problem with [13]. The scheme uses 
the Crank-Nicolson scheme for the time variable, and the spatial variable is 
discretized by the tension spline scheme on a non-uniform Shishkin mesh. 

Recently, in [19] Crank-Nicolson method for the time derivative on the 
uniform mesh and the central difference scheme for the spatial derivative on 
the Shishkin-type meshes is used for the problem in [13], and to enhance the 
order of convergence Richardson extrapolation technique is used. Also, in [7] 
singularly perturbed delay parabolic reaction-diffusion problem with mixed 
type boundary condition is solved. The problem is discretized by the implicit 
Euler method on a uniform mesh in the time and the extended cubic B-spline 
collocation method on a Shishkin mesh in the space variable. 

In real-world applications, higher-order methods are preferred to their 
lower-order counterparts since they provide better accuracy at a lower com- 
puting cost. This paper aims to design a uniformly convergent scheme with 
a higher order of convergence in both space and time variables or singularly 
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perturbed time delay reaction-diffusion problems. To discretize the problem, 
a similar approach to [14] is applied, that is, the Crank-Nicolson scheme for 
the discretization of the time derivative, and for the spatial variable, a hy- 
brid scheme on a special type of Shishkin mesh is used, which is gradually 
condensing starting from the center of the domain to both the right and left 
boundary layers. 

In this article, we considered the following singularly perturbed time- 
delay parabolic partial differential equation of the reaction-diffusion type of 
the form: 

O 


(a +L.,)u(x,t) = —b(a, thu(a,t — t) + f(x, t), (x,t) € D (1) 


subject to the initial and interval conditions: 
u(x,t) = dp(a,t), for (2,t) E To, (la) 
and boundary conditions: 
u(0,t) = dr(t) on T) = {(0,t) :0<t <T}, (1b) 


u(1,t) = ¢4,(t) on, = {(1,t):0<t<T}, (1c) 


where 
Leww(@,t) = —EUze (x, t) + a(x, t)u(z, t), 


0 < e << 1 and tT > 0 is the delay parameter. Also, D = OQ, x Q = 
(0,1) x (0, 7] and lf =T; UT, UT,, Ty = [0,1] x [-t, 0], Ty and T are the 
left and the right sides of the rectangular domain D corresponding to 7 = 0 
and x = 1, respectively. The functions a(z,t), b(x,t), f(a, t), dp, d, and ¢,. 
are assumed to be sufficiently smooth and bounded and satisfying, 


a(x,t) > a> 0,b(a,t) > B > 0, (x,t) € D. 


For the existence and uniqueness of the solution of (1), see [1]. Under the 
above assumptions the solution of (1) exhibits boundary layers along x = 0 
and «= 1. 


The article is organized as follows: In Section 2, we provide some prop- 
erties of the analytical solution and its derivatives. In Section 3, temporal 
semi-discretization, spatial discretization, and the derivation of the scheme 
are discussed. In Section 4, we investigate the uniform convergence of the 
fully discrete scheme. Numerical results and discussion are presented in Sec- 
tion 5. Finally, in Section 6, the conclusion of the paper is provided. 


Notation : Throughout this paper, C’ denotes a generic positive constant 


that is independent of ¢ and mesh sizes. Also ||.|| denotes the standard 
supremum norm, which is defined as ||f|| = Sup |f(a,t)|, for a function f 
x,t)ED 


defined on some domain D. 
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2 Analytical results 


In this section, the analytical aspects of the solution of problem (1) and its 
derivatives are studied. 


Lemma 1 (Maximum principle). Let ®(z,t) € C?+(D) and assume that 
®(z,t) >OonT. 


Then (2 4+ Lae) W(t) > 0,V(z,t) € D imply that ®(z,t) > 0,V(a,t) € D. 


Proof. Let the point («*,t*) € D, such that ®(x*,t*) = min ®(z,t) and 
assume that ®(2*,t*) < 0. Clearly, (x*,¢*) ¢ [ these implies that (a*,t*) € 
D. 

The function (x,t) attains its minimum at (x*,¢*), then 6, = 0,®, =0 
and ®,, > 0 at (a*,t*) . Therefore, from (1), it is easy to observe that 


0 
ay x ® it < ’ 
(2 + Le2)P(e,t) <0 
which is a contradiction as (4 + Le) (x,t) > 0. Hence (x,t) > 0 for all 
(x,t) € D. 


Lemma 2. The solution u of (1) satisfies the following bounds, 


Otsu 


Ox Oti |< Ce-V?,  for0<i+2j <8. (2) 


Proof. See [13]. 


In the proof of the error analysis, sharper bounds of the solution and its 
derivatives are required so that the solution u is decomposed into a regular 
component v, and a singular component w, as follows: 


Uu=vuU+uU, 
for more detail see [1, 18]. 
Lemma 3. The solution u of (1) satisfies the following bounds, 


Oitsy i 
llama! SCO +e? *”), 
— || < Ce~*? (exp(-Vax/Vé) + exp(—Va(1 — x)/ve)), 


for (a,t)€ D,O<14+2j7 <8. 


Proof. The proof of this Lemma is described in [1]. 
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3 Numerical scheme 


3.1 Temporal semi-discretization 


On the time domain [0,T], we use uniform mesh with step size At, such 
that" = tig = JAR Are 2 torg =O, lag MT Say t= mi, 
where y and m, are positive integer, and M is the number of mesh intervals in 
[0, 7]. To descretize the time variable for problem (1), we use Crank-Nicolson 


method, given by 


uI(a) = $o(x,—t;), for j =0,1,2,....m1,0<2 <1, 
(1 + Lon Jue) = (:- 46.0) wla) + At [P(e + FI*1(z)/, 


uit1(0) = @, wt!) =¢,, for j =1,2,...,M —1, 

(5) 
where F(x) = f(a,t;) — b(a,t;)u(z,tj-m,) and, wt" (x) = u(x, t;+41) is the 
semi-discrete approximation to the exact solution u(x,t) of (1) at the time 
level t; = jAt. The local truncation error of the semi-discrete method (5) is 
given by 

ej41 = u(a,tj41) — (a), 
where &(a) is the solution evaluated after one step of the semi-discrete 
scheme (5) taking the exact value u(z,t;) instead of wu? as the initial data. 
Re-write (5) in the form 


I+ Lon) we) = gi(x),2 € Qn, 


wt*(0)= gi, Wt(1) =o, 


(6) 


where g/(x) = 4! (Pw + FIt l(a) — L-nti(e)) + u(x,t;). 


For the stability of the Crank-Nicolson method, one can refer to [6, 8]. 


u(z,t))<C, (2,t)€D, i=0,1,2,3, the 


Lemma 4. Suppose that | 5p 


local error associated to scheme (5) satisfies: 
llejtill < C1(At)®, 9 =1,2,...,M. (7) 
Proof. Using Taylor’s series expansion, 


u(x, tj41) = u(x, t;) 


a = uz(z, t341/2) + O((At)?), 


= —Lent(a,tj41/2) + F(a, tj41/2) + O((At)*), 


(8) 


where F(a, t;41/2) = —b(a, t341/2)u(2, t341/2-m,) + f(2,tj41/2) and 


Iran. j. numer. anal. optim., Vol. 13, No. 2, 2023,pp 262-284 


267 


Fitted scheme for singularly perturbed time ... 


f(x, tj41) + F(z, t;) 


F(@,ty41/2) = + O((At)’), 


2 
Ga ba alt) | oaty?y, 
b(a, tj 41/2) = Metis) eae) O((At)?). 


From the above, we obtain 


u(x, tj41) ar u(x, t;) 
2 


Lew, tj41/2) = Lew + O((At)?). 


Clearly, the local error ||e;+1|| is the solution of the following BVP: 


(1 + Aan Jose = O((At)), 
€;+1(0) = 0, €;41(1) = 0. 


(9) 


Next, using the maximum principle for the operator (1 + 410.0) proves 


the result, for more detail please refer [6]. 


All the results in the above insure the second order convergence of the 
scheme (5). The global error of the time semi-discretization is given by 


J 
E; = u(az,t;) — w(x) = Ser. 
k=1 
Theorem 1 (Global error estimate). The global error estimate at t; is 


[les || A, fe 1D oy A 


Proof. see [11]. 


3.2 Spatial discretization 


In this section, we discretized problem (5), using a hybrid scheme. First, we 
define a special type of Shishkin mesh to discretize the domain 2, then the 
required scheme can be constructed. 
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3.2.1 Shishkin mesh 


For N > 2°, for integer ki > 2, the fitted mesh a. is constructed by dividing 
the domain Q,, = [0,1] in to three subintervals such that 


Q, = [0,0] U[o,1—o] U[1— a, 1]. 


Consequently, the mesh points and the transition parameter o are defined as 
described in [14, 26]. Let, L = L(N) satisfying, 


In(dnN)<L<inN and exp(—L)< 
Then, the location of the transition point is 
_ fl 
a= ming Znvel}, (10) 


where 7 > 4/,/a. Moreover, mesh points in the spatial direction is given by 


a for i= 1,2,..., %, 
m=§ otr(y—1/4)?+40(% 1/4), for P+17+2,.49 (1) 
1—ay-i, fori =241,242,...,N. 


where the coefficient r is determined from xy /2 = 1/2 . 

Let h max = max hi, where hy = 2; — 41,1 = 1,2,...,N, from (11) it 
L 

is clear that the maximum mesh width hma, always correspond to h x and 


Anyy of the domain Qy, that is, 


R max = hw =h fet < CN7}. 


wl2 


In the coarse part [7,1 —o], the mesh width satisfies the following, 


hina ON, (12) 
|hiza — hal < CN-?, (13) 


for more detail see [26]. 


3.2.2 Derivation of the scheme 


We use a hybrid scheme which is a combination of the central difference 
scheme and the fourth-order compact difference scheme to approximate the 
semi-discrete problem (5), where the coefficients are determined to make the 
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scheme exact for polynomials up to degree four and satisfy the normalization 
condition, gy +qf+q/ =1, fori=1,2,...,N—-1. 


ryjtl — .— Fit 1 j+1 i 
LNUP =r UP tr2 07 +97 ORY =O@) =a 9 ted tag. 


(14) 
where 


and the coefficients 7; ig? re ere FF dy , gf and at will be determined later. The 


development of the method is based on computing the local truncation error 
as follows: 
Tigins = ON Et) — £8 (OF) 


= LN Get) - a((r -£ Sean) 


15 
_ pn pith pe it atl — g-(qiti 4 Abe git oe 
=, (U1 TU + rf jig — G (aj) + Ey e,cUj_1) 
At At 
~j+1 ~jt+l1 ~jt+l1 ~jtl1 
> gf (ust + ria ee y= qi ( a a Lewitt a), 


Let’s denote the local truncation error corresponds to the discretization of 


the spatial variable by T*,, and T; zi+1 represent, Tj qit1 = AtTe.. 


Tea = Tpit Ta +... + TPa? + ORinax)s (16) 
where 
T? = (75, + Ti4 +75) — (qa = ar ety qi a a 
T} = (hig? — hiF,;) — (higsaf al fy — hegy af*)), 
mk aad 1 = 1 

T? = (hihi, Tr her 5) —2e- (2,,ata aN + hig qi a) 

k k-2 k 
ok — Mitts + ( ene gr ( HH Pisa gitl 
oe TU Rta TH \ aye 7 Ry Oe 


The truncation error is said to be of order p if Ti; = O(hijax) aS hmax — 0, 
for i= 1,2,...,N —1. The method is constructed “idee the conditions, 


Iran. j. numer. anal. optim., Vol. 13, No. 2, 2023,pp 262-284 


Amsalu Ayele, Andargie Tiruneh and Adamu Derese 270 


Tr =0, k=0,1,2, 


17 
iy = OGG a) k= 3,4. ( 


Therefore, the local truncation error JT’; can be written in the form, 


. mn ke 1 Za ee E Se\e9 
a Tu Tr Tu" 7 aaa > her, J+ ay Mn = hq; ) w 


1 os oe E SA Neg 
+ (Fins ~ Agr )+ qa + hig; a. 


(18) 


As a result, to fulfill the condition stated in (17), the coefficients are deter- 
mined as follows: 


~~ —2E — jt 
Mig ~ (a + hiss) + % M1), 


a+ + j+1 

r,, = | ——— +e &_1 |, 

as Cra tha)" ') (19) 
2 ; 

= (Ge + ata"), 1<i<N-1. 
aft 


The coefficients ar q; and q; ,i =1,2,...,N — 1 are defined in two different 
ways, 

(i) For the inner layer region, i.e (0,7) U (1 — a,1), the coefficients are given 
by 


ae hes 
tC Giga ea) 
2 
qi == _— (20) 
6 6Riai(y + Pigs) 
bes h? + h2,, + 38hihias 
ti Ghihixi 


(ii) For the outer layer region [0,1 — oc], depending on the maximum step 
length(hmaz) and ¢€, the coefficients are defined in two different cases. Define 
a=(at+ x); is a positive constant independent of ¢ and At, 

Case-1 When yh?2,,x||al| < ¢, the coefficients g;*,qf and q;, are defined by 
(20) 

Case-2 When yh2,,,.||G|| >, the coefficients g*,¢qf and g;, are given by 


max 


G =), =k, = 0. (21) 
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Remark 1. The coefficients qj, gf and q;_ are determined by the fourth order 
compact difference scheme in the boundary layer region, and in the regular 


region when yh?2,,,,|\@|| < ¢. While in the regular region these coefficients 


are determined by central difference scheme when yh?2,,,.||a|| > ¢. For more 
detail see [3, 14, 8, 5, 24, 25). 


3.2.3 Stability and error analysis 


To estimate the «— uniform convergence of the proposed scheme, we decom- 
pose the approximate solution of (14) into regular and singular components, 


U=V+W. 
Now, the truncation error of the method (14) can be decomposed as 
[Zi aoa | S Tepes | [ZG got | (22) 


where |T; 5i+1| and |T; ji+1| are the errors corresponding to the regular com- 
ponent 07+! and the singular component w/*!, respectively. 
Lemma 5. Let N > No, where No is the smallest positive integer such that 


An? Né 
=F Mllall + 2/At) < (23) 


*2? 
where L* = L(No). Then for any N > No, the coefficients 
Ti < 9, rg oO; ae <0O and regtrigtri; >0, for 7=1,2,...,N-—1. 


Proof. From (20) and (21) it is clear that g7 >0, gf >0, and gj > 0, then 


to show rj; <0, ri <0 and rf; > 0, we have 


_ At —2e = jt+l 4 2 
= tq, | a — 
"ii 2 be hes) Gi i-1 TAG) 
At —2e ; 2 
+ = bart gett 24 
"ig = 9 —< ay ee (off 7 x)): ve) 


a At Qe é jt1 2 
"og C= ae (« : At} )}" 


Using (20) and (21), the mesh width and the condition (23) proves rj; < 0, 
oe <0 and rf; > 0, then from (24), we obtain 


robe +r, = 0, for i= 1,2,...,.N—1. 


Remark 2. From Lemma 5 it is clear that the matrix associated with the 
discrete operator LY defined in (14) is an irreducible M-matrix and so has 
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a positive inverse. The operator LY in (14) satisfies the following maxi- 
mum principle on DN’. Hence, the scheme (14) is uniformly stable in the 
maximum norm. 


Lemma 6. Let ® satisfies 6 > 0 on PN. Then LNO >0on DN™ implies 
that ® > 0 at each point of DN™. 


Proof. The proof follows using Lemma 5 and Remark 2. 


Lemma 7. Let %"+! be the solution of (5) and U"*? be the solution of 
discrete scheme (14). Then the global error satisfies 


[art (x,) — UP *4|| < CAt(L/N)*. (25) 


Proof. For uniform mesh i.e o = 1/4, classical analysis can be used to prove 
the convergence of the scheme. So, we only consider the case 0 = o0,/éL. 
Case-1 Inner region [0,0] U [1 — 0,1], hi = hiti = ooVEN~1L, from the 
truncation error in (18), we have 


Pyrtt 


dx® 


iu 


Ten < cate( [is = hi|(hiqt + hi)? 


[vi—1,2441] (26) 


darth 
he) || 
+ ( a + hist) dx6 ) 
[vi—1,2441] 
Using Lemma 3, we can get 
Tis < CAt(L/N)* and Tig < CAt(L/N)* (27) 


Therefore, (22) gives 
Teg < CAt(L/N)*. 


Case-2 Outer region [c, 1 — a], determined by the relation between hmax and 
E 
(i) When yh?2,,x||@|| < 


max 


Part 
ies cAte( Iti — alata + ha)? dx® 
: xv 
[zs—1,2441] 
r ij dart 8 
+ RE + Pi) ae] - 
[ei-1,0541] 


Using (3), (12) and (13), we obtain the bound of the truncation error with 
respect to the regular component v is, 


Te < CAIN. (29) 


1,0 > 


Similarly for the singular component w using (4), (12) and (13), we obtain 
the bound of the truncation error with respect to singular component w, 
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since 
(exp(—Vax/Vé) + exp(-Va(1 — 2)/Vé)) < Cexp(-a(z, t)L), 


and 
exp(—L) < L/N. 


(ii) When yh2,,,|a|| > e, Assume L~* < CAt, 
Bartt 


dx? 


arti 
h2+h2 eae 
+( Fae 41) drt 


[wi—1,0441] 


Ta < Cte (Ihinihi 


a) 


(31) 
Using Lemma 3, (12) and (13), we obtain the bound of the truncation error 
with respect to regular component v and layer component w as follows 


Tis < CeAtN~?. 


Since hoa, > €; 
Tes < CN~* (32) 


and 


d? 
T®. < CeAt||—” 
ia 5 da? 


< CAt(exp(—Vax/Ve) + exp(—Va(1 — x)/Ve)) 
(ei-1>"i41] (33) 
< CAt(L/N)*. 
Hence, from (32), (33) and (22), we get 
Tig < CeAt(L/N)*. (34) 


Therefore, taking the estimate in case (i) and case (ii) then the truncation 
error estimate of the scheme (14) is given by 


Ti; < CeAt(L/N)*. (35) 


Finally, combining the uniform stability result and (35), concludes the proof. 


4 The full discrete scheme 


A combination of time and the spatial semi-discretization gives the full dis- 
crete scheme. Let U? be the numerical approximation for u(x;,t;) then the 
fully discrete scheme on the mesh D™ jis of the form 
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U4 (a;) = b»(2i, -t;), for j = 0,1,2,. ott and O51, 2 case N 


eel al OC ae pyle + 6 UI + rt - uit i =9, Gl, +¢G+g Gi, 

Us" = = = (t+ ‘i oon = “dott j+1)s 

for i=1,2,..,N—1 andj =1,2,...,M—1. 
(36) 

Where 
ij 
LN yitig _N U Ui =z Fi Xi 4 Fit zi), 
e074 6,04 At 
- Atl. 
Gf = Uf + | Pies) + Fa) - Ui, 


and the coefficients r; ;,17f;, Ta q; .4; and a. are described in subsection 
B22: 
The next theorem is for the ¢— uniform convergence of the proposed scheme 
(36). 


Theorem 2. Let u be the exact solution of (1) and U? be the approximate 
solution of (36). Under the hypothesis of Lemma (5). Then the global «— 
uniform error estimate of the scheme (36) satisfies the following bound 


||(ai, tj) — U2 || gv.m < o( (ar)? + (£/N)'), 1<i<Nand1<j<M, 


where L is as described in section 3.2.1. 


Proof. Case-1 For t < t, the delay term u(z,t — T) is not dependent on «, 
since the right hand side of problem (1) is —b(x,t)¢(a,t — 7) + f(a, t). Let 


a = OY x OM, where Q™ denotes the uniform mesh elements of [0,7], 


—N,M 
then the global error estimate on D_’ can be computed as follows 


||u(wi,ty) — US ll gyi S |lu(wi, ts) — Wl game + lid — G3] gum a 
+ ||O? - 


N,M. 
Uf llpe 


The combination of (7) and Lemma 3.2.3 result, 
\|u(a,, ¢ 7 U? Ihe < \|U2 — Uf |lppyam +car( (as)? + (L/N)*). (38) 
Taking the stability of the fully discrete method, we conclude that 
||? — Uf llggar S |u(wi, ty-1) — UP" gear (39) 


Now from (38) and (39) results a recurrence relation for the global error, so 
that 
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|lu(ai,t;)-U || 5" < c((ar?+(2/N)'), 1<i<Nand1<j <M. (40) 


Case-2 When t >, in this case the delay term u(x,t —) depends on ¢. The 
proof is done following the approach given in [1, 25]. The solution U(2;,t;) 
found on a is denoted by U;(2x;,t;).To determine the estimate of the 


=N,M 
global error on Dj, = QN x O38, where 034 denotes the uniform mesh 


with m mesh elements used on [t, 27], 

aaa (it fe U?)| = = 4; ht. Oe ee 1,tj- me) i (itp eae) +q, E; 
+ Sb! (Unt jr Whee) = u(vi, tien) + gq; E; 

Tq; ae (i tj-m,) — u(%i41, ieee) a qi E; 

+ Ti, 


(41) 


where E; and T; 4%; are the local error associated to the Crank-Nicolson 
method and the method in (14) respectively. Then using the bound (40), 


ie’ (u(ws,t;) — U2)| < ||Byll + [Teael + llee(westy) — Vell 
< o((ar? de win), 
1<i<Nand1<j<M. 


By introducing barrier functions and applying the discrete maximum princi- 
—N,M 
ple over Dj, , we found 


||u(ai,tj;) — Uillpyn < C( (At)? + (L/N)4), 1<i<Nand1<j<M, 


Lastly, by applying an induction argument we can get the required estimate. 


5 Numerical results and discussion 


To test the performance of the proposed scheme, we applied it to three differ- 
ent problems of the form (1). The maximum absolute errors (point-wise) and 
the rate of the convergence are calculated and tabulated. For those problems 
whose exact solution is unknown, the maximum absolute error is calculated 
using the formula: 


N,M _ N,M 2N,2M 
E =e ae (4,t;) -T0 (x2i, ta, )|- 
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where UN“ (x;,t;) denote the numerical solution obtained by N mesh inter- 
vals in the spatial direction and M mesh intervals in the time direction, such 
that M = T/At. The corresponding rate of convergence rate is given by the 


formula: 
N,M on at 
Re if = log2 pene - 
€ 


EN.M 


Also, the e— uniform maximum point -wise error is given by the for- 


mula: 
ENM — maxEN™ | 
€ 


and the corresponding ¢—uniform rate of convergence RN is given by the 


formula: 
N,M ee 
R™” = log2 ppaNam }° 


Example 1. Consider the following reaction-diffusion problem [9]: 
1 
Ut — EUgg + —u(z,t) = —2e-tu(e,t—1)+ f(z,t), (@,¢) € (0,1) x (0, 2], 


2 
moore, Gaeta, 


u(0,t) =e" + Pa aa) u(1,t) = e7 (+2) +e, te 0,2]. 


E (G2) 
The exact solution of problem (1) is u(x,t) =e (+) +e (e+ ve ), 


Table 1: Computed E™:4t and RN:4* for Example 1 


el N 64 128 256 512 1024 
At > 0.25 0.25/2 0.25/27 0.25/23 0.25/24 

2° 7.4762e-04 —1.5843e-5 3.7666e-5 —9.3469e-06 —2.3224e-06 
2.2385 2.0725 2.0107 2.0027 

2-2 1.0284e-03 2.4731e-04 —6.1734e-5 1.5411e-5 3.8514e-6 
2.0560 2.0022 2.0021 2.0005 

2-4 7.8737¢-4 1.9537e-4 4.8763¢e-5 1.2199e-5 3.0495e-6 
2.0108 2.0024 1.9990 2.0001 

2-8 — 5.6837e-4 1.4103e-4 3.52360-5 8.8022¢-6 2.2004e-6 
2.0108 2.0009 2.0011 2.0001 

2-8 —5.6713e-4 1.4042e-4 3.5088¢-5 8.7639¢e-6 2.1905e-6 
2.0139 2.0007 2.0013 2.0003 

2-10 5.8228e-04 1.4152e-04 3.5095e-05 8.7677e-06 2.190 7e-06 
2.0407 2.0117 2.0010 2.0008 


2-29 5 .8228e-04 1.4152e-04 3.5095e-05 8.7677e-06 —2.1907e-06 


2.0407 2.0117 2.0010 2.0008 
EN: 5.8228¢e-04 1.4152e-04 3.5095e-05 8.7677e-06 2.1907e-06 
RN: 2.0407 2.0117 2.0010 2.0008 
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Table 2: Computed EN:4* and RN»4t for Example 1, to show the order of convergence 
of space variable 


el N 32 64 128 256 512 
At > 0.25 0.25/4 0.25 /4? 0.25/43 0.25/44 

20 7.4762e-04 3.7666e-05 2.3324e-06 1.4570e-07 9.1115e-09 
4.3110 4.0134 4.0007 3.9992 

a-4 1.0284e-03 6.1735e05 3.85 14e-06 2.407 1e-07 1.5047e-08 
4.0582 4.0026 4.0000 3.9998 

oO 7.8752e-04 4.8772e-05 3.0501e-06 1.9063e-07 1.1915e-08 
4.0132 3.9991 4.0000 3.9999 

g-8 5.7001e-04 3.5282e-05 2.2070e-06 1.3793e-07 8.6217e-09 
4.0140 3.9988 4.0001 3.9999 

o-8 5.9289e-04 3.6752e-05 2.2962e-06 1.4374e-07 8.9838e-09 
4.0119 4.0005 3.9977 4.0000 

Q-10 1.0075e-03 6.2749e-05 3.9321e-06 2.4595e-07 1.7804e-08 
4.0050 3.9962 3.9989 3.7881 

g-12 1.1750e-03 7.2942e-05 4.5593e-06 2.8753e-07 1.7959e-08 
4.0098 3.9999 3.9870 4.0009 

2-14 1.1750e-03 7.2942e-05 4.5593e-06 2.8753e-07 1.7959e-08 
4.0098 3.9999 3.9870 4.0009 

2=16 1.1750e-03 7.2942e-05 4.5593e-06 2.8753e-07 1.7959e-08 
4.0098 3.9999 3.9870 4.0009 

EN:t 1.1750e-03 7.2942e-05 4.5593e-06 2.8753e-07 1.7959e-08 
RN:t 4.0098 3.9962 3.9870 3.7881 


u(x) 


(a) 6 = 2-8, N = 64, At = 0.25 (b) ¢ = 27-16, N = 64, At = 0.25 


Figure 1: Surface plot of numerical solution of Example 1 
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Figure 2: Log-log plot of Example 


) 


Example 2. Consider the following reaction-diffusion problem [1, 12, 13]: 


1+27 
2 


Ut — EUge + u(x,t) = —u(xz,t-—1)+t°, (2,t) € (0,1) x (0,], 
with 


u(x,t) = 0,(a#,t) € [0,1] x [-1,0], u(0,¢) =0,u(1,t)=0, te (0, 2}. 


Table 3: Computed EN“ and RN:4* for Example 2 


Method «|  N—128 256 512 1024 2048 
t > 0.25 0.25/2 0.25/22 0.25/23 0.25/24 
present 10° 2.7387e-04 6.7253e-05  1.6808e-05 4.2019e-06 _1.0505e-06 
2.0258 20004 2.0001 2.0000 
10-2 1.6900e-02 4.2067e-03  1.0448e-03 2.6493e-04 —9.2519e-05 
2.0062 2.0095 1.9795 1.5178 
10-4 1.9293e-02 4.8267e-03  1.2067e-03 3.0157e-04 —9.3848e-05 
1.9989 2.0000 2.0005 1.6841 
10-6 1.9357e-02 4.8428e-03  1.2067e-03 3.0275e-04 —9.3848e-05 
1.9989 2.0048 1.9999 1.6895 
10-8 1.9358e-02 4.8432e-03  1.2111e-03 3.0280e-04 —9.3833e-05 
1.9989 1.9998 1.9999 1.6902 
10-10 1,9358e-02 4.8432e-03  1.2111e-03 3.0280e-04 —9.3833e-05 
1.9989 1.9998 1.9999 1.6902 
10-12 1,9358e-02 4.8432e-03  1.2111e-03 3.0280e-04 —9.3833e-05 
1.9989 1.9998 1.9999 1.6902 
EN.t 1,9358e-02 4.8432e-03 1.2111e-03 3.0280c-04 9.3833e-05 
RN.t 1.9989 1.9998 1.9999 1.6902 


[12] ENE 2.39e-01 1.28e-01 6.60e-02 3.35e-02 1.71e-02 


[13] EN.t 2.78e-01 1.42e-03 7.13e-02 3.58e-02 1.79e-02 


Example 3. Consider the following reaction-diffusion problem [13, 12]: 
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Table 4: Computed E™:4t and RN:4* for Example 2, to show the order of convergence 
of space variable 


Pa N 332 64 128 256 512 
t > 0.25 0.25/4 0.25/42 0.25/48 0.25/44 

10° —-3.3960e-04 —2.1010e-05 1.3131e-06 8.2068e-08 —_5.1299e-09 
4.0147 4.0001 4.0000 3.9998 

10-2 .2.1214c-02 1.3267c-03  8.2918c-05 —5.1824c-06 —_3.2390c-07 
3.9991 4.0000 4.0000 4.0000 

10-4 2.4127e-02 —1.5092e-03 -9.4330e-05 5.8957e-06 —3.6848e-07 
3.9988 3.9999 4.0000 4.0000 

10-& 2.4192e-02 1.5136e-03 -9.4608e-05 5.913le-06 — 3.6957e-07 
3.9984 3.9999 4.0000 4.0000 

10-8 2.4192e-02 1.5137e-03 9.4615e-05 —5.9136e-06 —3.6960e-07 
3.9984 3.9999 4.0000 4.0000 

10-10 2.4189e-02 —-1.5137e-03  9.4615e-05 5.9135e-06 —3.6970e-07 
3.9982 3.9999 4.0000 3.9996 

EN:t 2.4192e-02 1.5137e-03 9.4615e-05 5.9136e-06 3.6970c-07 

RN.t 3.9982 3.9999 4.0000 3.9996 


(a) ¢ = 1077, N = 64, At = 0.25 (b) e = 10-4, N = 64, At = 0.25 


Figure 3: Surface plot for numerical solution for Example 2 
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Figure 4: Log-log plot of Example 2 
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Uy — EUze = —2eu(a,t—1), (x,t) € (0,1) x (0, 2], 


with 
u(0,t)=e"*, u(1,t) =e" VV), t € (0,2), (42) 
u(a,t) =e OVV®), (2, t) € [0,1] x [-1,0]. 
u(x,t) = e~+*/V®) is the exact solution of the problem. 
Table 5: Computed EN:At and RN:4* for Example 3 
Method eJ N +128 256 512 1024 2048 
At — 0.25 0.25/2 0.25/27 0.25/28 0.25/24 
present 10° 3.9168e-05 8.2891e-05 1.9790e-05 4.9090e-06 = 1.2254e-06 
2.2404 2.0665 2.0112 2.0022 
10-2 6.7565e-04 1.6837e-04 4.2083e-05 1.0517e-05  2.6290e-06 
2.0047 2.0003 2.0005 2.0001 
10-4 6.8750e-04 1.6903e-04 4.2110e-05 1.0518e-05  2.6291e-06 
2.0241 2.0051 2.0013 2.0002 
10-® 6.8750e-04 1.6903e-04 4.2110e-05 1.0518e-05  2.6291e-06 
2.0241 2.0051 2.0013 2.0002 
10-8 6.8750e-04 1.6903e-04 4.2110e-05 1.0518¢-05  2.6291e-06 
2.0241 2.0051 2.0013 2.0002 
10-19 6.8750e-04 1.6903e-04 4.2110e-05 1.0518e-05 2.6291e-06 
2.0241 2.0051 2.0013 2.0002 
10-1? 6.8750e-04 1.6903e-04 4.2110e-05 1.0518¢-05  2.6291e-06 
2.0241 2.0051 2.0013 2.0002 
EN:t 6.8750e-04 1.6903e-04 4.2110e-05 1.0518e-05 2.6291e-06 
RN: + 2.0241 2.0051 2.0013 2.0002 


[12] ENt 9.29e-03 4.43e-03 2.17e-03 1.07e-03 5.32e-04 


[13] EN 1.49e-02 7.75e-03 3.95e-03 2.00e-03 1.00e-03 


The fitted mesh technique applied in these computations, described in 
subsection (3.2.1), is uniform in the inner region, whereas in the outer region 
the mesh gradually condenses more and more from the center to both left 
and right ends of the interval. In Example 3 because of the boundary values 
there is no boundary layer on T’,., even though we use the same mesh tech- 
nique for this computation, it is possible to reshuffle the mesh appropriate 
to each specific example, for further description of this special characteristic, 
see [1]. In Tables 1 and 2, 3 and 5, for examples 1,2 and 3 respectively, 
the maximum absolute error and rate of convergence for different values of 
perturbation parameter ¢ and mesh numbers are presented, in all the case 
we observed that the computed e— uniform errors EN“! decreases monoton- 
ically as the number of mesh points increases, this ensures that the proposed 
scheme converges ¢— uniformly. The orders of convergence RN:“! is also in- 
dependent of the perturbation parameter ¢. To visualize the numerical order 
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(a) € = 1077, N = 64, At = 0.25 (b) ¢ = 10-4, N = 64, At = 0.25 


Figure 5: Solution graph for Example 3 
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Figure 6: Log-log plot of Example 3 


of convergence, the maximum point-wise errors are plotted in log-log scale in 
figures 2, 4 and 6. The effect of the perturbation parameter on the boundary 
layer behavior of the solution is also shown in figures 1, 3 and 5. The numer- 
ical result displayed in tables 2 and 4 reflects the actual theoretical order of 
convergence of the spatial variability of the proposed scheme. 


6 Conclusion 


Singularly perturbed time delay parabolic partial differential equation of 
reaction-diffusion type of the form (1) is considered. To get an approximate 
solution to this problem, we used a hybrid scheme which is a combination 
of a fourth-order compact difference scheme and a central difference scheme 
based on a special type of Shishkin mesh in space variable, and the Crank- 
Nicolson method on uniform mesh in the time variable. The proposed scheme 
is «— uniform convergent of order O((At)? + (L/N)*). Further, numerical 
experiments are carried out to demonstrate the performance of the proposed 
scheme. The numerical results clearly show the high accuracy and order of 
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convergence of the proposed scheme as compared to results available in the 
literature. 
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